function [b, b_Err, ycurve] = twoBodyDecayFit(xdata0,ydata0,ydata0_Err, averageFlag,xcurve0)

if max(ydata0) > 1e6             % fit density decay
    xdata = xdata0*1e-3; %[s]
    ydata = 1./ydata0*1e11;
    xcurve = xcurve0*1e-3;
    % Define fit function
    fun =@(b,x) b(2).*x+1/b(1);  
    % Define default initial guess
    b0 = [1.3, 1.5];   
    if averageFlag  % Fit with weight
        ydata_Err = (ydata0_Err*1e-11).*ydata.^2;
        W = ydata_Err.^(-2);
        [b,residual,jacobian,~,~,~] = nlinfit(xdata,ydata,fun,b0,'Weights',W);
    else
        [b,residual,jacobian,~,~,~] = nlinfit(xdata,ydata,fun,b0);
    end   
    % Calculate 95% confidence interval of the fitted parameters
    ci = nlparci(b,residual,'Jacobian',jacobian);   
    b = b';
    b_Err = b - ci(:,1);
    ycurve = 1e11./fun(b, xcurve);
    b(1) = b(1)*1e11;
    b(2) = b(2)*1e-11;
    b_Err(1) = b_Err(1)*1e11;
    b_Err(2) = b_Err(2)*1e-11;
else             % fit Number decay
    xdata = xdata0*1e-3; %[s]
    ydata = 1./ydata0;
    xcurve = xcurve0*1e-3;
    % Define fit function
    fun =@(b,x) b(2).*x+1/b(1);  
    % Define default initial guess
    b0 = [max(ydata0), 1/max(ydata0)/0.6];   %[ N(t=0), 
    if averageFlag  % Fit with weight
        ydata_Err = (ydata0_Err).*ydata.^2;
        W = ydata_Err.^(-2);
        [b,residual,jacobian,~,~,~] = nlinfit(xdata,ydata,fun,b0,'Weights',W);
    else
        [b,residual,jacobian,~,~,~] = nlinfit(xdata,ydata,fun,b0);
    end   
    % Calculate 95% confidence interval of the fitted parameters
    ci = nlparci(b,residual,'Jacobian',jacobian);   
    b = b';
    b_Err = b - ci(:,1);
    ycurve = 1./fun(b, xcurve);
    b(1) = b(1);
    b(2) = b(2);
    b_Err(1) = b_Err(1);
    b_Err(2) = b_Err(2);   
end

end